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ABSTRACT Dark matter, believed to be present in many galaxies, is interpreted as a 
hydrodynamical system in interaction with the gravitational field and nothing else. An 
equation of state determines the mass distribution and the associated gravitational field. 
Conversely, the gravitational field can be inferred from observation of orbital velocities 
of stars in the Milky Way, in a first approximation in which the field is mainly due to the 
distribution of dark matter. In this approximation, the equation of state is determined 
by the gravitational field via the equations of motion. 

The potential is an exact solution of the equations of motion in the approximation 
of weak fields, 

a / , df 

-A</>ocp, 0oc— , 
dp 

where / is the free energy density. The second equation (the integrated hydrostatic 
condition) determines p in terms of <f>; the first equation then becomes a nonlinear 
equation of the Emden type that issolved exactly by the chosen potential. 

The resulting equation of state is a simple expression that accounts for the main 
features of the galactic rotation curve over 6 orders of magnitude. 

I. Introduction 

One of the enduring problems of astrophysics is to place an upper limit on the 
mass of a star. Let us agree from the outset that the observed mass of a spherically 
symmetric object is the asymptotic value of the function M that appears in the quasi- 
Schwartzschild metric, 

„ 2M(r)G. , x 

Soo W = c 2 l ^- . 1.1 

r 

A locally observed mass is defined in a region where this function is slowly varying. The 
Great Attractor near the center of the Milky Way is observed, at a distance from the 
center of around 10 16 cm, to have a value for this parameter that is several million solar 
masses (Ghez 2008). This is as much as 5 orders of magnitude greater than "reasonable" 
physical models (Hartle 1978) . 

Analysis of the distribution of velocities of orbiting stars show that the newtonian 
potential cannot be attributed to visible sources; the locally observed mass increases 
far too rapidly with the distance from the center. Both problems can be qualitatively 
explained in terms of 'dark matter', the high value of M because the equation of state 
of dark matter is unknown and not subject to the physical constraints of known forms 
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of matter, the unexpected variation of M with distance because dark matter may be 
present in regions that appear to be empty. 

All that is known about dark matter is that it does not interact with ordinary 
matter. It does not interact with electrodynamics, it is not in thermal equilibrium 
with ordinary matter or with radiation, and the temperature is not defined. This 
puts the theoretician in the same position as he confronts in hydrodynamics when the 
temperature is eliminated from the theory by means of the ideal gas equation and the 
equation of state reduces to a relation between density and pressure. The free energy 
density is a function of density alone, the entropy density df/dT is zero and the pressure 
is 

P=(P^-1)/(P). (1-2) 

The system is thus determined by the expression chosen for the function /(p); this 
expression, or the inferred relation between pressure and density, will be referred to as 
the equation of state. 

In this paper we shall propose a simple, analytic expression for the newtonian 
potential that accounts for the main features of the rotation curves of our Galaxy. 
From this we shall determine the unique equation of state that is required in order 
that Einstein's equations admit this idealized potential in the weak field approximation. 
Then we use the equation of state so determined in the full system of Einstein's equations 
in the presence of dark matter. 

If this equation of state turns out to be applicable in other galaxies as well, then 
this approach to the problem of dark matter can be considered as an alternative to 
modified gravity; see for example Delbourgo (2008) and Mannheim (2011). 

Some data. The radius of the Milky Way is about ro = 10 23 cm and the mass is 
about 2MG = 2 x 10 17 cm. The innermost, observed satellite has a nearest approach of 
about 2 x 10 15 cm and it moves in the newtonian field of a mass of about 2MG = 10 12 . 
The cgs system is used throughout; lkpc = 3 x 10 21 cm and 8nG = 1.863 x 10~ 27 cm/ g. 

Summary 

Our model for (the negative of) the gravitational potential is 

2MG . r ~~\~ b 52 
= Mr) = k In , b = e cm. 

r r 

The equations of motion include the hydrostatic condition in integrated form, 

2 r alp' 

where / is the free energy function. Einstein's equations, in the weak field approxima- 
tion, give a unique equation of state represented parameterically as follows 

f(p) = Bip sinh 4 ip — p, p = A sinh 4 -0, p = B J sinh 4 ip dip, 

A and B constants. We solve the relativistic equations of motion using this equation of 
state to obtain the gravitational metric and the density distribution of the Galaxy. 
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II. The equations of motion 

We shall calculate static, spherically symmetric solutions of Einstein's equations, 

G»„ = ^T„„, G = .7414xl0- 28 — , 
c z g 

with a metric of the form 

ds 2 = e^cdt) 2 - e A dr 2 - r 2 dO, g 00 = cV«, g rr = -e A «, 

and a matter energy momentum tensor of the form 

Too = pU U , T rr = pg 00 , (2.1) 

all other components zero. Besides Einstein's equation we invoke the hydrostatic con- 
dition in integrated form 1 

c 2 <9f 

¥'~-»=k- (2 ' 2) 

The reduced form of Einstein's equations given in the textbooks, beginning with 
that of Tolman (1934), is 1 



G/ = — e _A (— + -2] + ~2 = —SttGp/c 2 . < 2 ' 3 > 



With the notation 2 



K(r)=e^ = l-^. 

they are 



H(r) = e" A = 1 ^, K(r) = e" 



H' = ^— - - 8nGr P (e- u - p/c 2 p), 
K' = 8nGH- 2 rp. 

Some computer programs do not like very large numbers and work better if we 
change variables, introducing x by 

r = e x , 



1 See Section VI. 

2 The choice of the letter m in the first expression is traditional, but unfortunate, 
in as much the locally observed mass defined in (1.1) is 2M(r)G = m(r) + u(r). (See 
below.) 



to get 

—m{x) = w r 3 p — p/pc 2 ), 

dx ' , 2A s 

To continue we need an expression for the free energy that will allow us to express 
the density and the pressure in terms of the fields, with the help of Eq.(2.2). To deter- 
mine the free energy we shall work, provisionally, with the weak field approximation. 



III. The equation of state 

A weak field approximation will be used to determine an approximate equation of 
state, subject to later adjustment. In this approximation we replace Eq.s (2.4) by 

m'{r) = r 2 wp (3-1) 

K'(r) = rwp. (3.2) 

The primes, as before, denote the derivative with respect to r. The newtonian potential 
is —0/2, where 

1 - e"" « (1 - H) + (1 - K) = - + -=: 0. 

r r 

Combining (3.1-2) we get 

2 jj t i 2 

m = —r <p , wp = m /r 
The full set of equations is thus 

y^=^, -r- 2 (r 2 d>y = wp. (3.3) 

Something is known about 0, from observation of radial acceleration of orbiting 
stars. A family of satellites moving in circular orbits with radius r in a radial, newtonian 
potential V have orbital speed v given by v 2 = rV'. Observation has revealed that there 
is a wide interval in which the speed is nearly constant, independent of the distance, 
which implies that, in this interval, the potential is approximated by V = —0/2 = 
(k/2) ln(r). We shall model the function 4>(r), then calculate the equation of state. In 
other words, when the distribution 0(r) is known from observation, then the last pair of 
equations provides a parametric representation of the relation between the free energy 
density / and the density p. Finally we shall use this equation of state in the exact, 
relativistic field equations. 
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Example 

Taking 

= Mn^^, r = e x , (3.4) 
r 

we obtain the velocity distribution shown in Fig.l with k = 1, b = e 52 . It is very nearly 
constant for x < 52 and very nearly newtonian for x » 52. 

m = —r 2 (p' = kb (= kb h r > b), 

r + b r 



and the density 



wp(r) = —r (r cp ) 



2 , 2 , / \ / k 



r 2 (l + r/b) 



2 ' 
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Fig.l. The orbital velocity distribution that was used as a model of the observations 
in the Milky way. The velocity is shown in km/sec. 



The parametric representation of the equation of state is thus, in this case, 



2 9/ j , , , t + b k 
= 0(t) = fcln — — , wp[t) 



c 2 dp t ' rw t 2 (l + t/b) 2 ' 

Equivalently, 

16 k 

p= TT -sinh 4 ^, ii:=(j)/2k. (3.5) 
b z w 

The free energy is obtained by integrating the hydrostatic equation, 

A#=0 l^L = (j) ^p = 2k(—(^p)- 



Thus 

k2 



Ap) = i6 c 2fc2 -ffc) = y sinh4 tpdip = ip sinh 4 ?/> - p. 



The last term, 



p :-- 



is the pressure, 



sinh 4 ipdip = — ( sinh(4-0) - 8 sinh(2'0) + 120 
32 



16 k 2 c 2 A d , „ 
cr w op 



(3.6) 



(3.7) 



It is shown function of ip in Fig. 2. 
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Fig. 2. Equation of state. Upper curve: the normalized pressure p as a function of 
the variable ip, see Eq.(3.6). The lower curve is the n = 4 polytrope, p = (1/5) sinh 5 t(j, 
a perfect fit at low values of ijj. The innermost observed orbiter, at x = 35, is at if) = 7.5 
and the highest value (at x = 21.6) is 10.5. The deviation from the ordinary polytrope 
is considerable. 



Returning to the equations (3.1-2) - the weak field approximation - we now apply 
the equation of state in the form (3.5-6), fix the appropriate initial values, 



x = a = 52, m(a) 



bk 



.2 x 10 17 , K(a) = 1 -k In 2 + 



k 



kln2 



2 ' 2' 

and run the equations in Mathematica. The program runs from x = a in both directions, 
covering 40 orders of magnitude of the radius, correctly reproducing the exact solution 
(3.4). 



Remark. The function 



ip = - ln(l + b/r) 
2 



is an exact solution of the modified Emden equation 

Aip + sinh 4 -0 = 0, A = r~ 2 ^-r 2 ^-. 
b 1 dr dr 

The original Emden equation has ip n instead of sinh 4 ip; it has an exact solution in the 
case that n = 5 only. 
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IV. Solutions of the relativistic equations of motion 

Using the same equation of state we now solve the exact Einstein equations numeri- 
cally. Once the equation of state has been found there are no free parameters. 

Upgrading the equations, from the weak field approximation (3.3) to the exact equa- 
tions of motion (2.9) has limited effect, for the fields are relatively weak everywhere; 
that is, 4> « 1. 

After adjustment of the boundary values we obtained a solution covering the range 

16 < x < 77, 8.9 x 10 6 < r < 2.8 x 10 33 . 

The values of m(r) at some chosen values of the local mass, are 

m = 10 n , a: = 39.1, r 

m = 10 12 , £ = 41.4, r 

m = 10 13 , x = 43.7, r 

m = .4 x 10 17 , x = 58, r 



= 9.6 x 10 iO , 
= 9.5 x 10 17 , 
= 9.5 x 10 18 , 
= 1.5 x 10 25 . 



The zero of the function m(r) is a computational error. It was verified that the 
calculation gives result of high accuracy for 25 < x < 60. However, it is not m(r) 
that should be interpreted as the local mass, but r<p(r), since (p rather than m/r is the 
newtonian potential. For rcj) the corresponding values are 



rcj) = 


3 x 10 5 , 


X 


= 23.0, 


r = 


= 9.7 x 10 9 


rcj) = 


10 11 , 


X 


= 36.4, 


r = 


6.4 x 10 15 , 


rcj) = 


10 12 , 


X 


= 38.8, 


r = 


7.1 x 10 16 


rcj) = 


10 13 , 


X 


= 41.3, 


r = 


8.6 x 10 17 


rcj) = 


.4 x 10 17 , 


X 


= 58.0, 


r = 


.5 x 10 25 . 



Observation of the innermost satellites of the Milky Way suggests a local mass of 
about 10 12 (3 million solar masses) at a distance of 10 16 from the center. 

Other numerical results are as follows. The density is positive; there is a characteristic 
bump in the density profile - see Fig. 3b - where the density reaches the highest value, 
p = 100g/cm 3 at x = 21.6, 

Pmax = 100^/cm 3 at r = 2.4 x 10 9 . 

This "object" is comparable to our Sun, in size, mass and gravitational field strength. 
At the shortest distance observed for an orbiting satellite, r = 2 x 10 15 (x = 35.23), 
the density is about 1.3 x 10 -10 . The pressure has a similar profile - Fig. 4, with a peak 
value of 2.5 x 10 16 . The gravitational field — (p(x) = c~ 2 goo — 1 also has a maximum 
-see Fig. 5, reaching a maximum value of ep = .00003 at the same point. This is about 
ten times stronger than the gravitational potential at the surface of the Sun. The 
appearance of such shapes is very common when polytropic equations of state are used; 
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they are relativistic features not seen in the weak field approximation. For some stars 
the maximum value of </> can rise to get very close to the limiting value of unity, at 
which point a horizon would appear. 

The local mass predicted by the model at the distance of the inner orbiters is less 
than what is observed, by about one order of magnitude. If this discrepancy can be 
removed by refinements of the model, then we will have a picture of the galactic center 
that is very different from a Schwartzschild black hole. 




Fig. 3. A plot of In pj lnr against x = lnr. 
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Fig. 5. The pressure profile in the inner region.. 
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Fig. 6. The potential has a maximum at the same point. The density and pressure 
peaks are narrower because of the high value of the "polytropic index" . 
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Fig. 7. The velocity distribution predicted by the relativistic equations of motion, 
in km/sec. Compare Fig.l. 
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V. The nature of dark matter 

The equation of state was obtained in parametric form, Eq.s (3.5) and (3.6), 



kc 2 

p = Csinh 4 ^, p = C(sinh4V> -8sinh2^ + 12^), (5.1) 

with C = 16k/b 2 u. It bears a remarkable similarity to an equation first proposed by 
Stoner(1932) and used by Chandrasekhar (1935), 3 

p = Ci sinh 3 1, p = C 2 ( sinh At - 8 sinh 2t + 12t) . 

The basis for this formula is a model of fermions in a collapsed state, the Femi sea being 
filled up to q/m = sinh^. 

The similarity, if not regarded as a coincidence, suggests that dark matter may be 
a cloud of "ice crystals" that consist of fermions in a highly reduced state. The total 
absence of interactions, and of photons, is a premise of Chandrasekhar's work. A slightly 
different model reproduces our equations (5.1) exactly, but the fermions need to have 
an additional degree of freedom; for example, an extra dimension of momentum space. 
With q = sinh?/', E = cosh?/', integrating over the 3-sphere with radius q : 4 

4 1 f j^d A q ^ f 4j i 1 f _ d 4 q 



The equation of state that has been developed here is consistent with the observed 
velocity distribution, even for the innermost orbiters, but this does not give enough 
information to develop a microscopic model of dark matter. 

It is essentially a hydrodynamical system. At very low densities the star is an n = 4 
polytrope, 

p=-/5 5 / 4 , p = sinh 4 if;. 



With 
and 

it works out to 



5 

16k _ 24 
p = ap, a = —5 — = 6.1 x 10 , 

p = af3p, p = c 2 k = 9.4 x 10 14 , 



p = Ap b /\ A = -a" 1 / 4 /? = 6.0 x 10 20 . 
5 



3 It was 'corrected' and used by Oppenheimer and Volkov (1939) in their study of 
neutron stars. The original version, quoted here, is in Landau and Lifshitz (1958) page 
168. 

4 The factor E in the first integral arises because the mass of the compound system 
is the sum of the energies of the constituents. This factor cancels the factor 1/E in the 
volume element, just as in the calculations of Stoner and Chandrasekhar. 
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At very high densities the adiabatic index is effectively infinite, p = p/A and 

p = = 2.35 x 10 14 p- 

This relation is confirmed at the point where ift = 10.48, p = .00122(7 / cm, where we are 
in the regime of high densities. 

VI. Theory 

The reader will have noticed that the equations that were used are classical, but a 
closer look reveals some novelties. 

1. The integrated hydrostatic condition is classical; a short calculation shows that 
taking the gradients of both sides leads to the usual hydrostatic condition, 

p grad 4> = — grad p. 

The difference is that the integrated form incorporates the boundary condition that 
fixes the speed of light at infinity. 

2. It is of interest to ask why this boundary condition has not been applied previously. 
One part of the answer is that the stars of Eddington and Chandrasekhar all have abrupt 
boundaries at a point where the temperature and the density vanish, so that these 
fields are not continuous. Another part of the reason is that the equations employed by 
Eddington (1926), Chandrasekhar (1935) and many others differ in one particular from 
ours: The factor e~ u that multiples the density in Eq.(2.8) is absent; consequently, the 
function v is represented only by its derivative, so that fixing its boundary value has no 
meaning. 

In view of this difference we need to justify our approach. All equations used are 
variational equations based on the following action, 

A= 8^gJ d * x ^ R + J d A x^(p{g^^, v -c 2 )-f{pj). 

The non relativistic approximation of the matter lagrangian is one that was used by 
Fetter and Walecka (1980) to obtain a variational formulation of hydrodynamics: the 
equation of continuity and the Bernoulli equation (in integrated form). 5 The term f(p) 
is the free energy, as is seen from the structure of the equations of motion (Fr0nsdal 
2007, 2008) . The variational approach to hydrodynamics is an application of the Gibbs 
variational principle to the case that the temperature is frozen, so that neither tem- 
perature nor entropy plays any role, as is appropriate for a treatment of dark matter, 
effectively a hydrodynamical system. The density p is denoted p + p by Eddington; 
this is a matter of notation, and irrelevant in the immediate context, since p/c 2 p never 

5 The non relativistic approximation is taken by setting \I> = c 2 t + neglecting 1/c 2 
terms and interpreting $ as the velocity potential. 
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exceeds 10~ 4 . The gradient of the field ^> corresponds to Tolman's vector field U ; for 
a stationary solution \1/ o is a constant, while Tolman's normalization condition leads 
to Uq = ^goo- This is what gives rise to the cancellation of this metric function in 
Tolman's equations of motion, and it constitutes an important difference in principle 
between our approach and that of Tolman. 6 

3. Another consequence of action principle dynamics is that the current is conserved, 

The usual approach does not admit a conserved current and breaks with non relativistic 
theory in this respect. 

Although the current is conserved, it is not directly related to the "mass" as we would 
define it. Our approach to stellar structure is to start the analysis from the outside, 
using observational data. In the models considered here the metric has the asymptotic 
form 

_ 2 2MG 

c g 00 = l , 

r 

with M constant. This is an obervational datum, measured by observing the motion 
of test bodies. In this paper, that is what we call the mass of the star. Now it is 
always pointed out that the function m that appears in (2.4) tends to M at infinity. 
In the traditional approach the equation is just dm/dr = wr 2 p and the mass can be 
expressed as an integral over the density, provided that m vanishes at the origin. 7 
Our theory preserves the continuity equation of classical hydrodynamics, but there is 
no direct connection between mass and the conserved quantity. 8 Consequently, there 
is no need to postulate that the function m vanishes at the origin. 

An equation of state in hydrodynamics is a relation between density and pressure. 
The pressure term in Einstein's equation is not as important as the role that is played by 
the pressure in the hydrostatic equation. In our approach the equation of state follows 
from the expression that is chosen for the free energy density. Our approach preserves 
all the structure of hydrodynamics, including the equation of continuity. 

Ultimately, the model of our Galaxy must be improved by including the contribution 
of visible matter; the overall structure can be studied in terms of an idealized, contin- 
uous, spherically symmetric distribution that makes an important, additional contri- 
bution to the action. The absence of any interaction between dark and visible matter 

6 In a variational approach it is important to specify the independent variables; any 
constraint is the source of great complications. 

7 It has been pointed out that the integral J pd s r does not have the correct measure. 
Kippenhahn and Weigert (1990) calls this situation hazardous, but no one seems to 
have taken the warning seriously. 

8 It is true that the increment of mass, between distances r± and ri from the center, is 
the integral of the density over the region bounded by two spheres. Note, nevertheless, 
that the non relativistic gravitational potential arises entirely from the time compo- 
nent of the metric, while the function m is in the space component. The unfortunate 
association of this function with the mass is due to Tolman's normalization condition. 
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makes this very straightforward; in the absence of any interaction between the two kinds 
of matter the free energy density is additive. Note that it is essential to recognize the 
roles of two quite different density fields. 

VII. Speculating about the center 

In this paper the study of our galaxy has been approached from the outer regions, 
because that is where observations have been made up to this time. The analysis 
is especially interesting because nothing is known about the nature of dark matter. 
Consequently, there is no need to ask what pressures and densities can be allowed; 
there is no way that we can answer questions of this kind. 

Historically, a number of statements have been made that would place limits on 
the mass of certain types of stars. These statements all rely on two assumptions: 
(a) that the nature of the matter within the star is within the limits of our knowledge 
and (b) that the total mass is the integral of the density, from the center outwards. If 
these premises may said to be reasonable as far as stars made up of ordinary matter is 
concerned, they cannot be applied with any degree of confidence to the new situation 
that is faced in connection with dark matter. The bump in the density is a common 
feature of relativistic models; here it is predicted to occur at a distance of about 2.4 x 
10 9 cm from the center. It is conceivable that future observation may validate this 
prediction, but what happens inside is beyond our reach. 
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